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Abstract 

We have studied the phase diagram of two capacitively coupled Josephson junction arrays with 
charging energy, Ec, and Josephson coupUng energy, Ej. Our results are obtained using a path 
integral Quantum Monte Carlo algorithm. The parameter that quantifies the quantum fluctuations 

E 

in the i-th array is defined by ai = Depending on the value of Oj, each independent array 

may be in the semiclassical or in the quantum regime: We find that thermal fluctuations are 
important when a < 1.5 and the quantum fluctuations dominate when 2.0 < a. Vortices are the 
dominant excitations in the semiclassical limit, while in the quantum regime the charge excitations 
are important. We have extensively studied the interplay between vortex and charge dominated 
individual array phases. The phase diagrams for each array as a function of temperature and 
inter-layer capacitance where determined from results for their helicity modulus, T(a), and the 
inverse dielectric constant, e~^{a). The two arrays are coupled via the capacitance Cinter at each 
site of the lattices. When one of the arrays is in the quantum regime and the other one is in the 
semi-classical limit, T{T,a) decreases with T, while £~^{T,a) increases. This behavior is due to 
a duality relation between the two arrays: e. g. a manifestation of the gauge invariant capacitive 
interaction between vortices in the semiclassical array and charges in the quantum array. We find a 
reentrant transition in T(r, a), at low temperatures, when one of the arrays is in the semiclassical 
limit (i.e. ai = 0.5) and the quantum array has 2.0 < a2 < 2.5, for the values considered for the 
interlayer capacitance of Ci„ter = 0.26087,0.52174,0.78261,1.04348 and 1.30435. Similar results 
were obtained for larger values of a2 = 4.0 with Cinter = 1.04348 and 1.30435. For smaller values 
of Cinter the supcrconducting-normal transition was not present. In addition, when 3.0 < a2 < 4.0, 
and for all the inter-layer couplings considered above, a novel reentrant phase transition occurs 
in the charge degrees of freedom, i.e. there is a reentrant insulating-conducting transition at low 
temperatures. Finally, we obtain the corresponding phase diagrams that have some features that 
resemble those seen in experiment. 

PACS numbers: 74.81.Fa, 03.75.Lm, 73.43.Nq, 05.10.Ln 
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I. INTRODUCTION 



This paper considers the phase diagrams of two capacitively coupled Josephson junction 



arrays (JJA), made of ultrasmall 



many theoretica 



imental 




unctions. Two-dimensional JJA have been the subject of 



isl Il9 1 and exper- 







29 



30, 



14 



16 



17 



3l| studies. Advances in submicrometer 



technology |2y] and in nanolitographic techniques 



.es 



have allowed the fabrication of 



relatively large arrays made of ultrasmall superconductor-insulator-superconductor Joseph- 
son junctions (JJ). The J J areas may vary from a few microns to submicron dimensions, 
with self-capacitances ~ 3 x 10~^ fF and nearest neighbor mutual capacitance Cm ~ 1 
fF. Notice that the mutual capacitance can be at least two orders of magnitude larger than 
the self-capacitance. 

In these arrays there are two competing energies: the Josephson coupling energy, Ej, 

2 

and the charging energy Ec = of the junctions, due to the charge localization in the 

islands, where e is the electronic charge. In the semiclassical limit, Ej >> Ec, the phase 

of the superconducting order parameter of the junctions is well defined. The Josephson 

coupling induces fluctuations in the charge number that produces a supercurrent, where 

the average Cooper pair number is undefined. In this regime the vortex excitations are 

pinned by the intrinsic lattice potential and the array is in a superconducting state. In 

the quantum limit, Ej << Ec, the electrostatic energy to add one Cooper pair in one of 

the two islands in a junction is much larger than the thermal energy. The electric field 

localizes the Cooper pairs in the islands and the quantum fluctuations of the phase of the 

superconducting order parameter delocalizes the vortex excitations. This charge localization 

in the junction islands drives the array to an insulating state. The number of Cooper pairs, n, 

and the phase of the superconducting order parameter, 0, obey the Heisenberg's uncertainty 
1 

principle, AnA(f) > -. This uncertainty relation has been demonstrated experimentally in 
Al - AI2O3 - Al junctions |32 1 . 

The superconducting-insulating (S-I) phase transition, induced by the charging energy 
in arravs of this type, has been experimentally measured as a function of a by groups in 
Delft j3| and Harvard j^]. Their junction sizes had constant values, while they varied 
the normal state junction resistance to change the Josephson coupling energy. This allowed 
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them to fabricate arrays with a in the range [0.13 — 4.55] 23|, or as high as 33 j2J]. JJA 
have also been studied in connection with quantum phase transitions The JJA have the 
advantage, over films, that their internal structure can be carefully controlled experimentally. 
The drawback is that the array sizes are limited by fabrication constraints. There is a large 
literature studying the phase transition structure of one JJA with ultra small junctions as 
the temperature, the external magnetic field and a are varied (for a comprehensive recent 
review see 

A novel experimental system composed of two capacitively coupled JJA made of ultrasmall 
junctions was discussed .n a paper by the Delft group H. wlriclr each anay was produced 
with different a parameters. Initial theoretical studies of this system were started in j33|], [3J| 
and in j^^. In this paper we study the phase diagram of this very interesting system. The 
Hamiltonian describing the coupled arrays can be formally written as 7t! = Hj{l) +Hj{2) + 
i/c'(l,2), where the Hj{i) denote the Josephson Hamiltonians for each array, and i/c'(l,2) 
describes the two arrays capacitive interactions. Hc{^i 2) includes the total charging energy 
matrix, which includes the self- and mutual capacitive terms in each plane, plus the arrays 
interaction due to the ultrasmall nearest neighbor capacitive coupling between them. This 
Hanultouran . more complicated to =tudy than the single array problem studied before Q 
l25l |. Here each array is described by the ratio ai = (EcjEj.) [i = 1,2), with the mutual 
array capacitive coupling denoted by aint- In the a, -C 1 limit, the i — th array is dominated 
by localized vortex excitations, Vi, while the Cooper pair excess charge excitations, Qi, are 
in a superconducting state. In the aj ^ 1 regime, the system has the Qj's localized in an 
insulating state while the V^'s are delocalized. There are different parameter regimes that 
can be studied. Here we consider in detail the case when one array is V dominated and the 
other array Q dominated. Using the Villain approximation and duality transformations , 



it was shown l3c 



3J, |35[ that the interaction between vortices and charges has a minimal 



gauge coupling form, and that the interaction is sharply localized in space: A vortex and a 
Cooper pair interact only if they are located at the same site in both arrays. In previous 
studies an effective Hamiltonian was derived, in terms of Vi and Q2, which shows a high 
degree of dual symmetry, with the <^==4> Qo interaction expressed as a capacitively 
induced minimal gauge like coupling 

We used a path integral Quantum Monte Carlo algorithm to further study the behavior 
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of this system. We calculated the phase diagram described in terms of the vortex helicity 
modulus, T, and the charge inverse dielectric constant, e~^, as a function of a and tempera- 
ture. We find that as a varies there are different novel reentrant phase transitions seen in T 
and e^^. These re-entrances indicate transitions between SC-N, SC-I and I-N phases. The 
layout of the paper is as follows: In section |n] we define and explain the model of interest 
and the path integral formalism to evaluate its partition function, as well as the physical 
quantities that are calculated to characterize the behavior of the model. In section IIIII we 
present and discuss the results of our calculations and finally, in section IIVI we summarize 
our main findings and present our conclusions. 



II. THE MODEL AND ITS PATH INTEGRAL REPRESENTATION 

In this section we define the model that describes two JJA made of ultra small junctions 
capacitively coupled at each lattice site. The parameter values used in the simulations 
correspond to those present in the experimental samples fabricated and studied at Delft j^l ■ 
The dimensions of the experimental arrays were = 230 and Ly = 60, with intra-layer 
mutual capacitance of Cm = 2.3 fF, and inter-layer capacitance Cint = 0.6 fF. 

The model Hamiltonian is defined by j^^, 

^ = T ^ ^ Mn)C-^!M,r2)n,if,) + Ej, J2 - COS {Mn)- Mr 2))) 

<fi,r2> fJ-=^,'^=^ <ri,r2> 

+ Ej, (l-cos(02(ri)-02(r2))), (1) 

<ri,r2> 

where the sums are over nearest neighbors, and fi and u label the matrix elements of the 
mutual capacitive interactions, fl and denote the positions of the junctions in arrays 1 
and 2, respectively. The operator ri^(ri), for the number of Cooper pairs, and the phase of 
the superconducting order parameter 0jy(r2), satisfy the Heisenberg commutation relations, 
[ri^j(ri), 0,^(r2)] = —iSffj^^ff^^f^^u- The charge carried by the Cooper pairs is Q = 2e, and Ej-^ 
and Ej^ are the Josephson coupling constants within each array. C~^^^v is the electric field 
propagator and its inverse, C^^y, is the block matrix representing the geometric capacitance. 
In what follows we will use the notation C^,^ = C^^ to denote the diagonal capacitance 
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matrix. The block matrices are given by, 



—C 

^ Cint ) 

0, 



zC„,^^ + Cint) , if /i = and n = f2, 

if /i = z/ and r*i = r 2 ± m, 

if /i 7^ z/ and n = r2, 
otherwise. 



(2) 



Here d is a unit vector and z is the lattice coordination number. The off-diagonal blocks are 
written as —Cint I n,n with I the identity matrix and the size of the arrays. The block 
matrix elements in the Hamiltonian in Eq. are written as 



Ci,i - Ci ^ [I - CintCa ^] , 

C2,2 = C2 ^ [l - CintCi ^Cg ^] , 

C~ /~( — — 1 [t z^ — l/"! — 11 

1,2 — "-^int*-^! ^2 ~ '-"int'^1 ^2 J 

C" — 1/^ — 1 Ft /~i^l/^^ll 

2,1 — '-^mt*-^2 *-"l ~ '-"int'^2 *-"l J 



(3) 



To evaluate the partition function, S, we can calculate the trace over the phase operators 
or over the trace of the number operators fi. The path integral representation of H can be 



obtained using the states 



, .^.| exp[m(rl)0(r^)] 



^271 



(4) 



Following the steps outlined in reference [15|, we obtain the lattice path integral represen- 
tation of the partition function. 



L^-l 



27r 



27r 



exp< - 



/3 ft 



dr 



T=0 r n(T,r) 

+ i^n{T,f)^{r,r) + Ej ^ [l - cos (0(r, rl) - 0(r, rl)) 

f <fl,f2> 



(5) 



Eq. (0) involves the phase 0(r, r) and the charge integer fields n(r, r) as statistical vari- 
ables, in a three dimensional lattice with two spatial dimensions x Ly, and one imaginary 
time dimension, Lr. The angular phases (f){T, r) G [0, 27t] are defined at the nodes of the 
lattice with periodic boundary conditions in the x and y space dimensions. The integer 
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fields, n{T, r), lie in the bonds between two consecutive nodes along the imaginary time axis 
r, and they can take any integer values. The quantization condition in S is imposed in 
terms of the periodic boundary conditions 0(L^,r) = 0(0, r). The ^ oo limit has been 
formally taken to replace the sum over imaginary time by its integral. This H representation 
is amenable to numerical calculations in comparison to its operator representation. 

The arrays can be in a superconductor, insulator, or normal states, depending on the 
values of T, a and Cint- To characterize the superconducting or normal state behavior we 
calculated the helicity modulus, T, for each array. 



T 



9k2 



k=0 



(6) 



T measures the energy needed to carry out a phase twist between the boundaries of the 
array along the k direction. The helicity modulus is proportional to the superfiuid density 
per unit mass ps, 

,,.)4(!^)T(r). (T) 

Here a is the lattice spacing, m the Cooper pair mass, and V the volume. Combining Eqs. (jSj) 
and (ini), we obtain the path integral representation of the helicity modulus when the twist 
is along the x axis jisl ]. 

Lr-l 




r, fy) - 0(r, fl + x) ^ 



r, r^j - 01 r, + x 



r, r^ + x 




(8) 



T=0 fl 

The charge coherence in the arrays is determined by the inverse dielectric constant e~^, 
defined as. 



lim, 



1 



1 



-{n{q)n{-q)) 



Combining Eqs. (0) and 
charge number operator 
function (n{fi)n{f2)) as 

1 



El) , as well as Fourier transforming the capacitance matrix and the 
15| , we can write the path integral representation of the correlation 



n(r^)?T,(r^) 



■C(rl,rl) + f ^ C(rl,ri)C(r^,rl)/m(r5),m(rl) 



(10) 



Substituting this result in Eq.Q yields the inverse dielectric constant expression 



-C(A;)(|m(A;)| ) 



(11) 



In these equations the path integral representation of the integer fields m{r)'s is m(r*) = 
Ylt=o^ m(r, r). It is important to note that the path integral representation of in Eq. 
is not exactly the inverse dielectric function of a gas of Cooper pairs, since it depends on 
the discretization of the imaginary time axis in Lr slices. Nonetheless, we expect that it 
contains most of the relevant information of the dielectric properties of the gas of Cooper 
pairs in the arrays. 

III. RESULTS AND DISCUSSION 

A. Parameter values in the simulations 

To carry out the Monte Carlo moves in the phases we used the standard Metropolis 
algorithm. To speed up the calculations, we replace the U{1) continuous symmetry in the 
phases by a discrete subgroup Using = 5000 has proved to yield good results. 
Discretizing the phases has the advantage of using integer arithmetic that allows building 
cosine tables for the Boltzmann factors in the partition function. This simplification cannot 
be used for the integer fields, except when = 0, in which case the m(r, r) fields can 
be summed up exactly. This approach allows to build up look up tables that introduce an 



adequate effective potential 

Once the system reaches thermal equilibrium after Aether (between 10^ and 10^) MC 
sweeps, the thermodynamic averages of interest were calculated. A measure was taken after 
Nsweeps passes through the arrays updating the phases, and Mgweeps passes updating the 
integer fields. In the semi-classical limit, a <^ 1, and we needed to consider N sweeps = 1 
and Msweeps = 1 at high temperatures (T > 0.4), and N sleeps = 4 and Msweeps = 4 at low 
temperatures (T < 0.4). In the quantum limit, a > 1, we took Nsweeps = 4 and Msweeps = 4, 
at high temperatures, and Nsweeps = 10 and Msweeps = 10 at low temperatures. These 
parameter values were chosen to minimize the decorrelation times due to the long range 
charge interactions. Proceeding in this way we carried out averages over 16384 MC steps 
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at high temperatures and 32768 at low temperatures. Error bars were calculated using the 
biased reduction method described in reference 

We studied the helicity modulus and the inverse dielectric constant for each array as a 
function of temperature, a, and the ratio of the inter-layer self capacitance to the intralayer 

Cint 

mutual capacitance, -y^- The reduced temperature was varied in the range [0.05, 1.0] in 
0.05 steps at high T's and 0.025 steps at low T's. To study the transition between the 
semi-classical and quantum states, the quantum parameter in one array was kept fixed 
at «! = 0.5 while we varied the quantum parameter of the second array taking values, 
= 0.5, 1.0, 1.25, 1.50, 2.0, 2.5, 3.0, and 4.0. In addition, the capacitance ratios were integer 

Cint 

multiples of the Delft's experimental parameters 23|: — ^ = k x 0.26087, with k = 1,2,3,4, 
and 5. The values chosen for this ratio allowed us to study the effects of the capacitive 
coupling between the arrays, going from weak coupling, to the strong coupling regime when 
K > 3. The simulations were carried out in lattice sizes L^. x = 32 x 32 and Lj. = 32, 
for 0.5 < a < 2.0, and up to 96 for higher values 2.0 < a < 4.0. The larger the value a 
the larger the imaginary time axis had to be. For the Lr values chosen we obtained reliable 
results, with negligible finite size effects along the imaginary time axis. 



B. Results for T and e ^ 

In this section we present and discuss the results for Ti^2 and e^l for different values of 
the quantum parameter 0^2 while ai = 0.5 was kept fixed. We start considering the case 
when both arrays are in the semi-classical regime, ai = 0.5 and 0:2 = 1.0 and 0.26087 < 
Cint ^ 1.30435. The corresponding results are not shown explicitly here, however, we briefly 
describe them to compare, where appropriate, with previous studies. We found that at 
low temperature Ti 2 are finite -SC phase- with Ti > T2, with small thermal fluctuations. 
Both quantities decreased monotonically down to zero as temperature increased. For T > 
^sc(«i,2) they were equal to zero indicating that the arrays were in the normal phase. The 
superconductor to normal (SC-N) transition temperature, Tgc'(a), shifted downwards as 
the charging energy increased. Around this transition temperature the fluctuations of T 
became larger. Since the arrays were in the semiclassical regime the charges did not play 
an important role in the behavior of the system, except for the small decrease in phase 
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coherence, and we obtained, = 0, in the whole temperature range. These results for the 
two coupled arrays are in agreement with previous studies for one array The 
results described here indicate that Tsc{o:) was weakly affected by the interlayer capacitive 
coupling. 

In Fig.^we show the helicity modulus -upper panels- and the inverse dielectric constant 
-lower panels- of each array as a function of temperature when ai = 0.5 and 02 = 2.5, 
for low and high values of The dependence of the SC-N phase boundary Tscicn) on 
Cint will be shown and discussed later. In figure ^ we see that the second array starts to 
develop charge coherence: At very low temperatures the inverse dielectric constant of 1, 
although small, it has increased while the phase coherence of the second array decreased, 
as shown by the reduced value of the helicity modulus. This behavior confirms the dual 
behavior between vortices in one array and charges in the other one, due to their gauge 
interaction, as predicted in reference p&l\. There, an effective action was written in terms of 
four interacting imaginary time Coulomb-like gases. The effective Hamiltonian was shown 
to be dually symmetric between charge and vorticity with complicated kernels. In the 
simplest case, where one array has one vortex and the other one charge, their interaction 
has a minimal gauge coupling that is proportional to the site coupling capacitance between 
the arrays. As a consequence of this gauge capacitive interaction, the fluctuations in the 
helicity modulus and the inverse dielectric constant are much stronger than those observed 
in the semi-classical limit in a single array. Nonetheless, the quantum fluctuations in the 
second array are not sufficiently strong to significantly affect the semi-classical behavior of 
the first array. This indicates that the gauge interaction between vortex and charge Coulomb 
gases is weak. On the other hand, it also found that as increases the charge coherence 
of the second array decreases, restoring its phase coherence. In addition there appears to 
be a reentrant phase transition in the second array at low temperatures, as suggested by 
the concave curvature of its helicity modulus and the slight convex curvature of its inverse 
dielectric constant. This picture becomes more evident when increasing further the charging 
energy in the second array, to ^2 = 3.0, while keeping ai = 0.5. The results are shown in 
Fig. 121 When = 0.52174, as the system cools down, array 2 undergoes a SC-N phase 
transition at T ^ 0.35, followed by a reentrant N-SC transition at T ~ 0.25. In addition, 
we also observe that at T ~ 0.45 the charge coherence increases significantly and then 
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decreases at T ~ 0.15. This suggests an insulating to normal reentrant transition. Once 
again, this is a manifestation of the dual character between the charge and phase degrees of 
freedom in the two arrays as well as the gauge interaction between the vortex and charge 
Coulomb gases. As the interaction capacitance increases further to ^ = 0.78261, the 
reentrant transition in the phase degrees of freedom is less evident, but it is still present. 
A similar situation is observed in the charges degrees of freedom. Increasing further the 
arrays interaction capacitance improves the phase coherence of the arrays at the expense of 
their charge coherence. Note that the thermal fluctuations are strong in the phase and in the 
charge degrees of freedom for the array parameters considered here. In spite of the capacitive 
interaction between arrays 1 and 2, the behavior in array 2 does not affect significantly the 
behavior of array 1. When the quantum parameter of array 2 is sufficiently large, 0:2 = 4.0, 
the quantum effects become dominant, and the charge degrees of freedom destroy the phase 
coherence when = 0.52174 and 0.78261. This behavior is shown in figure 01 For this 
value of a2 we needed to consider Lr = 96 for array 2 to get meaningful results. We found no 
transition at all in the phases, but we do get a sharp N-I transition for the charges, although 
the fiuctuations in those degrees of freedom are rather large. Increasing Cjnt further leads 
to a different scenario. For = 1.04348, a SC phase appears in the temperature range 
0.25 <T< 0.45, followed by an insulating phase at lower temperatures. For = 1.30435 
we obtain a similar situation. In this case the phase coherence is greater than for the previous 
values of the interaction capacitance parameter. These results are shown in figure 01 

As pointed out in sectional when the quantum parameter, a > 1, the numerical eval- 
uation of the path integral representation of e^^ -Eq. ()11|) - depends importantly on the 
time slice discretization of the imaginary time axis, that is, there are finite size effects. To 
ascertain this situation we calculated T(T) and e~^(T) for four different imaginary time axis 
lengths, Lr = 48,64,96, and 128. In figures a), (b) we show the results of these calcu- 
lations for Tj(T), z = 1,2 when = 1.30435 with ai = 0.5, and 02 = 4. We see that 
Ti(T) shows no finite size effects since the data for different values of L^- for each T fall on 
top of each other. There are very small finite size effects for those temperatures close to 
the SC-N transition temperature. However, T2(T) shows important finite size effects in the 
temperature region where the reentrant transition sets in. This reentrance becomes better 
defined, and its behavior becomes smoother for Lj. > 96. Similarly, the SC-I temperature 
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shifts slightly to lower temperatures for > 96, while the N-SC temperature remains about 
the same. On the other hand, ej~^(T) = in the whole temperature range considered. This 
is not unexpected since array 1 is in the semi-classical regime {ai = 0.5), and the role of 
the charge degrees of freedom is negligible. Because of this we do not show e^^^ in figures |S1 
Nonetheless, e2^(T) is sensitive to the Lr length axis, as can be seen in figure E^c). In this 
case, we see that the insulating phase sets in at the temperature boundary Tsc-i{Lt) ■ For 
Lr > 96 it shifts slightly to lower temperatures. We also see strong fluctuations in the in- 
verse dielectric constant at about Tqc-i- For other values of and 02, the results are less 
significant. Although we did not see a clear saturation of 63 ^(T) at the reentrance transition 
for the two largest values of Lj. we were able to minimize the finite size effects along the 
imaginary time length. Carrying out calculations using larger imaginary time axis becomes 
very demanding from the computational point of view. 



C. Phase diagrams 



Here we present the cumulative results from the estimation of the critical temperatures 
from our T(T, a) and e~^(T, a) calculations that provide the phase boundaries for the SC-N 
and I-N phase transitions. In array 1 the quantum parameter ai = 0.5, is kept fixed, while 
a2 was varied in the interval 1.0 < 02 < 4.0, in Aa2 = 0.5 steps, for each one of the values 
of Cint/Cm- The transition temperatures were estimated from the behavior of T and 
as functions of temperature, for each value of 02- In the classical arrays, according to the 
Berezinskii-Kosterlitz-Thouless (BKT) theory, the critical SC-N temperature occurs at the 
intercept of the helicity modulus with a straight line with slope |. In the quantum arrays 
it has been shown experimentally that in the limit Cg/Cm —>■ there is a crossover from a 



conducting to an insulating phase 



2i 



29[. Due to the finite size of the experimental 



arrays there is a rounding of the transitions, because the screening length is shorter than the 
size of the arrays [2^. Theoretically it has been argued that for any finite screening length 
the phase transition is washed out even in arrays of infinite size i39\. This happens because 
in the BKT scenario the SC-N transition crucially depends on the unscreened nature of 
the logarithmic interaction between vortex pairs j36|,|40|. Nonetheless, it has been recently 
shown jl8, 3| that even in the regime of strong quantum fluctuations the data for T(T), at 
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low temperatures, can be very well fitted to the Kosterlitz's renormalization group equations 
in a finite size analysis of T. In spite of the theoretical and the experimental differences 
in estimating the transition temperatures, here we will use the BKT theory since we have 
not carried out a detailed finite size analysis. When there is a low temperature reentrant 
phase transition, the transition temperature is determined from the temperature at which 
the physical quantity of interest vanishes, as is done with the arrays's electric resistance 
in experiments. The error bars in the transition temperatures were estimated taking into 
account the size of the temperature variation step, AT, in the calculations. Since AT = 0.05, 
then the error bars in the transition temperatures of the phase boundaries are 5Tc = 0.05. 

In figure 10 we show the phase diagrams T^ versus 02 for the phase and charge degrees 
of freedom for arrays 1 and 2, when Ci^t/Cm =0.52174. We see that the SC-N temperature 
boundary of the semiclassical array 1 remains almost unchanged as a function of 02, with 
Tc ~ 0.85. The SC-N phase boundary of array 2, however, changes significantly as 02 
increases. The transition temperature decreases monotonically and it extrapolates to zero at 
a2 = 4. At this 02 value the second array is fully in the quantum regime. Similar results were 
obtained for Cint/Cm =0.26087. It was also found that the I-N phase boundaries of the arrays 
are qualitatively similar for 1 < 0:2 < 2.5 with Cint/C*m = 0.26087 and 0.52174. That is, for 
array 1 Tic (02) = 0, and for array 2 the I-N transition temperature increases monotonically 
and reaches a maximum at 02 = 2.5. However, when array 2 is in the quantum regime, 
a2 > 2.5, the situation is different. There we found that the I-N boundary of array 1 is equal 
to zero when Ci^t/Cm = 0.26087. Increasing this quantity up to 0.52174, yields a nonzero 
and monotonically increasing Tic(a2)- For array 2 with Ci^t/Cm =0.26087 we found that the 
^ic(«2) decreases slowly and for 3 < 0:2 < 4 it becomes constant, Tjc — 0.3. For larger values 
of the interaction parameter, Cint/C*m = 0.52174, the I-N phase boundary decreases down to 
zero for 3.0 < a2 < 4.0, indicating the existence of a reentrant novel phase transition. The 
results for Ci„t/C„ =0.52174 are plotted in figureEl For Cint/C^ = 0.78261, the SC-N phase 
diagram for array 2, and the I-N phase diagram of both arrays, are modified slightly. The 
Tsc-N boundary of array 2 decreases at about the same rate when 1.0 < a2 < 3.0 for the 
values of Cint/C*m given above. Nonetheless, the transition temperature for 0:2 = 4.0 is not 
zero, instead, it moves upwards. Tsc-n may be approximated by a straight line for the whole 
interval of 02 considered here. The I-N boundary, it increases monotonically from 02 = 1.25 
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up to = 2.0, where it reaches a maximum and then slowly decreases appearing to reach 
a minimum at 0:2 = 3.0. Above this value it increases slowly. These results are not shown 
here since they are similar to those in Fig. [7| A further increase in the coupling capacitance 
between the arrays to, Cint/C*m = 1.04348, does not change the SC-N and the I-N boundaries 
of array 1, but it changes the SC-N phase diagram for array 2, and the I-N phase diagram of 
both arrays. The SC-N transition temperature of array 2 decreases linearly as a function of 
0^2, while the I-N boundary increases monotonically in the temperature range 1.0 <T< 2.5. 
It reaches a maximum at T = 2.5 and then decreases slowly up to 0:2 = 3.0, where it has a 
minimum and then increases again up to 0:2 = 4.0. Note that the position of the maximum is 
shifted to higher values of a2, and the transition temperature at 02 = 4.0 shifts downwards 
compared to the previous value of Cint/C'm- When Cmt/Cm is at the highest value considered 
here, 1.30305, the SC-N boundaries of arrays 1-2 become straight lines. The former has zero 
slope and the latter shows a negative slope, i.e. the SC-N transition temperature of array 2 
decreases linearly as 02 increases. On the other hand, the I-N boundary of array 1 starts as 
a horizontal straight line for 1 < ^2 ^ 1-5 then increases reaching a maximum at ^2 = 2.5, 
then decreasing monotonically down to zero at 0:2 = 4.0. These results suggest, again, the 
existence of a reentrant transitions in array 2 when it is in the quantum regime. In contrast, 
the I-N boundary of array 2 is a horizontal line at T = 0, for 1 < 02 < 3.0, and Ti_N(a) 
becomes nonzero for 02 > 3.0. 

The structure of the phase diagrams indicate that by increasing the capacitive coupling 
between arrays and increasing the 02 parameter one can induce novel reentrant phase tran- 
sitions not only in the phase degrees of freedom but also in the charge freedoms. This is a 
situation that does not occur in single JJA. 

IV. CONCLUSIONS 

We have carried out extensive path integral quantum MC simulations of two capacitively 
coupled JJA with ultra small junctions. One in the semiclassical limit and the other in the 
quantum regime. We studied the behavior of the helicity modulus and the inverse dielec- 
tric constant as a function of temperature for different values of the interlayer capacitances. 
When both arrays are in the semiclassical regime, regardless of the interlayer coupling consid- 
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ered here, each array shows a SC-N transition at finite temperatures. For these semiclassical 
arrays the charge degrees of freedom contribution is neghgible and there is no I-N transition 
at finite temperatures. As array 2 enters the quantum regime (2.0 < 0:2 < 2.5), a SC-N 
reentrant phase transition appears and the fiuctuations in T2 and become significantly 
larger due to quantum fiuctuations as explained in references Qj- At the same time 
the quantum array starts to develop charge coherence with an finite temperature I-N tran- 
sition. This scenario occurred for all the values of the interlayer coupling considered here. 
As array 2 becomes more quantum (3.0 < 0L2 < 4.0), and for all the interlayer couplings 
considered, reentrant transitions occur not only in the phase degrees of freedom, but also in 
the charges degrees of freedom. When 02 = 4.0, the SC-N transition is washed out for val- 
ues of the interlayer coupling C\nt/Cja < 0.78261, having only an I-N transition. Increasing 
further the interlayer coupling yields a SC-N reentrant phase transition at low tempera- 
tures for 0.20 ^ T < 0.5, together with the usual I-N phase transition for the charges. 
These scenarios can be understood as a manifestation of the gauge interactions between the 
phase degrees of freedom and the charge degrees of freedom that are present in each array. 
The latter results depend on the values of the quantum parameters. Finally, we showed a 
series of interesting phase diagrams for both arrays where we found some resemblance to 
the reentrant-like behavior observed experimentally by the Delft's group in 2D JJA with 
charging energy [2^ [si | . 
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FIG. 1: T and e ^ as a function of temperature, ai = 0.5 (circles), for the semiclassical array, and 
a2 = 2.5 (triangles), for the quantum array. The values of are (a) 0.52174 and (b) 0.78261. 
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FIG. 3: Same as in Fig.dwith a2 = 4.0 and %^ = (a) 0.52174 and (b) 0.78261. 
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FIG. 5: Lt finite size behavior of Ti^2 and e2 for Lj- =48 (circles), Lr =64 (triangles), Lj- =96 
(squares), and Lt =128 (diamonds). The interaction between arrays is = 1.30435. (a) semi- 
classical array (ai = 0.5), (b) and (c) quantum array (a2 = 4.0). e^^ is equal to zero in the whole 
temperature range and is not shown. 



22 




FIG. 6: Estimated transition temperatures of each array versus ol^ for = 0.52174. Array 1 
(•), array 2 (A). The solid lines represent the SC-N phase boundary while the dashed lines the 
1-N phase boundary. 
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FIG. 8: Same as in Fig.Elwith ^ = 1.30345. 
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